Introduction to R and Python
Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο GitHub
1 Πρώτα βήματα
1.1 Εγκατάσταση βιβλιοθήκης
Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)
1.3 Εγκατάσταση βιβλιοθηκών
Με την εντολή packages('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν
Δοκιμάστε το
1.4 Φόρτωση βιβλιοθηκών
Εντολή για φόρτωση των πακέτων
1.6 Αλλαγή του working directory
Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]
1.7 Εισαγωγή αρχείου
Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε
Όπως αντιλαμβάνεστε, σήμερα θα δουλέψουμε με δεδομένα τα οποία ήδη υπάρχουν εντός της R. Όμως για την εργασία που θα παραδώσετε εσείς, θα χρειαστεί να φορτώσετε δεδομένα από ένα αρχείο excel, το οποίο θα σας δοθεί στο τέλος της άσκησης. Για να το κάνετε αυτό, θα χρειαστεί να χρησιμοποιήσετε την βιβλιοθηκη readxl και την εντολη read_excel1
Ας δούμε τα δεδομένα μας:
## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
| Sepal.Length | Sepal.Width | Petal.Length | Petal.Width | Species |
|---|---|---|---|---|
| 5.1 | 3.5 | 1.4 | 0.2 | setosa |
| 4.9 | 3.0 | 1.4 | 0.2 | setosa |
| 4.7 | 3.2 | 1.3 | 0.2 | setosa |
| 4.6 | 3.1 | 1.5 | 0.2 | setosa |
| 5.0 | 3.6 | 1.4 | 0.2 | setosa |
| 5.4 | 3.9 | 1.7 | 0.4 | setosa |
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
## [1] setosa setosa setosa setosa setosa setosa
## [7] setosa setosa setosa setosa setosa setosa
## [13] setosa setosa setosa setosa setosa setosa
## [19] setosa setosa setosa setosa setosa setosa
## [25] setosa setosa setosa setosa setosa setosa
## [31] setosa setosa setosa setosa setosa setosa
## [37] setosa setosa setosa setosa setosa setosa
## [43] setosa setosa setosa setosa setosa setosa
## [49] setosa setosa versicolor versicolor versicolor versicolor
## [55] versicolor versicolor versicolor versicolor versicolor versicolor
## [61] versicolor versicolor versicolor versicolor versicolor versicolor
## [67] versicolor versicolor versicolor versicolor versicolor versicolor
## [73] versicolor versicolor versicolor
## [ reached getOption("max.print") -- omitted 75 entries ]
## Levels: setosa versicolor virginica
Και ας μετονομάσουμε τα ονόματα των ειδών2:
1.8 Λιγότερα δεκαδικά
Με την ακόλουθη εντολή, θα εμφανίζονται πλέον λιγότερα δεκαδικά:
1.9 Στατιστικά στοιχεία για κάθε taxon
Θέλουμε να μάθουμε την ελάχιστη, την μέγιστη και την μέση τιμή, καθώς και την τυπική απόκλιση3:
## $`Iris setosa`
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## Min. :4.3 Min. :2.3 Min. :1.00 Min. :0.10 Iris setosa :50
## 1st Qu.:4.8 1st Qu.:3.2 1st Qu.:1.40 1st Qu.:0.20 Iris versicolor: 0
## Median :5.0 Median :3.4 Median :1.50 Median :0.20 Iris virginica : 0
## Mean :5.0 Mean :3.4 Mean :1.46 Mean :0.25
## 3rd Qu.:5.2 3rd Qu.:3.7 3rd Qu.:1.57 3rd Qu.:0.30
## Max. :5.8 Max. :4.4 Max. :1.90 Max. :0.60
##
## $`Iris versicolor`
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## Min. :4.9 Min. :2.0 Min. :3.0 Min. :1.00 Iris setosa : 0
## 1st Qu.:5.6 1st Qu.:2.5 1st Qu.:4.0 1st Qu.:1.20 Iris versicolor:50
## Median :5.9 Median :2.8 Median :4.3 Median :1.30 Iris virginica : 0
## Mean :5.9 Mean :2.8 Mean :4.3 Mean :1.33
## 3rd Qu.:6.3 3rd Qu.:3.0 3rd Qu.:4.6 3rd Qu.:1.50
## Max. :7.0 Max. :3.4 Max. :5.1 Max. :1.80
##
## $`Iris virginica`
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## Min. :4.9 Min. :2.2 Min. :4.5 Min. :1.40 Iris setosa : 0
## 1st Qu.:6.2 1st Qu.:2.8 1st Qu.:5.1 1st Qu.:1.80 Iris versicolor: 0
## Median :6.5 Median :3.0 Median :5.6 Median :2.00 Iris virginica :50
## Mean :6.6 Mean :3.0 Mean :5.6 Mean :2.03
## 3rd Qu.:6.9 3rd Qu.:3.2 3rd Qu.:5.9 3rd Qu.:2.30
## Max. :7.9 Max. :3.8 Max. :6.9 Max. :2.50
## Iris setosa Iris versicolor Iris virginica
## Sepal.Length 0.35 0.52 0.64
## Sepal.Width 0.38 0.31 0.32
## Petal.Length 0.17 0.47 0.55
## Petal.Width 0.11 0.20 0.27
Μπορείτε επίσης να χρησιμοποιήσετε το σύνολο εντολών από την πρώτη εργαστηριακή άσκηση, προκειμένου να δημιουργήσετε τα στατιστικά αυτά δεδομένα. Δοκιμάστε το για το IQR (Interquantile Range).
1.10 Γραφική απεικόνιση των αποτελεσμάτων
Ας αναπαριστήσουμε γραφικά τα δεδομένα μας για κάθε taxon και μια μεταβλητή4:
## Violin plot
ggplot(iris, aes(Species, Petal.Length)) + geom_violin(adjust = 4) + geom_boxplot(width = 0.1,
fill = "black", outlier.colour = NA) + stat_summary(fun.y = median, geom = "point",
fill = "white", shape = 21, size = 2.5) + xlab("Taxon") + ylab("Petal Length") +
theme(axis.title = element_text(face = "bold"), axis.text.x = element_text(face = "italic",
colour = "black"))## [1] "G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/EKPA/2021"
Κάντε το ίδιο με τις υπόλοιπες μεταβλητές και αλλάξτε το όνομα του y άξονα.
ΠΡΟΣΟΧΗ
Για την εργασία σας, θα χρειαστεί να ελέγξετε εάν οι μεταβλητές σας έχουν κανονική κατανομή και εάν εμφανίζουν μεγάλο συντελεστή συσχέτισης. Ανατρέξτε σε προηγούμενες εργαστηριακές ασκήσεις για να δείτε πώς μπορείτε να το κάνετε αυτό.
2 Κυρίως ανάλυση
Είμαστε πλέον σε θέση να προχωρήσουμε στην κυρίως ανάλυση των μορφομετρικών μας δεδομένων. Είναι το πιο δύσκολο και κρίσιμο βήμα και θα χρειαστεί να είστε εξαιρετικά προσεκτικοί όταν θα τρέξετε την ακόλουθη εντολή για την εργασία σας:
attach(iris) ## For your assignment, change that to mydata
iris <- iris[apply(iris,1,function(row) all(!is.na(row))),] ## Do the same here
set.seed(2934) ## Change that if you want to whichever number you like
cv.grp <- sample(1:10, nrow(iris),
replace = TRUE) ## Remember to change the name
cv.lda <- c() ##A function that splits our data to 2 independent subsets, one to
## train our model and one to test how well our model behaves
for(i in 1:10) {
train.ind <- which(cv.grp!=i)
test.ind <- which(cv.grp==i)
model.train <- lda(Species ~ Sepal.Length +## Change the variable names
## according to the ones that
## appear in your data
Sepal.Width + Petal.Width + Petal.Length, ## Do the same here
data = iris, subset = train.ind) ## Change to mydata
cv.lda[test.ind] <- predict(model.train,
newdata = iris[test.ind,])$class ## Do the same here
}
mean(cv.lda==as.numeric(iris$Species)) ## And here## [1] 0.98
Lec.lda <- lda(Species ~ Sepal.Length + ## This is the test point
Sepal.Width + Petal.Width + Petal.Length, ## Change the variable names
data = iris) ## Change to mydata
plot(Lec.lda,col = c("green","blue","red", "black", "orange")
[iris$Species], pch = 19) ## Change to mydataΜε την ακόλουθη εντολή βλέπουμε τα αποτελέσματα της ανάλυσης μας:
Lec.lda
## Let's say what we did
##--------------------------------------------------------------------------------##
# Call: lda(Species ~ Sepal.Length + Sepal.Width + Petal.Width + Petal.Length,
# data = iris)
##--------------------------------------------------------------------------------##
## How much of the data lie in each species
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.33 0.33 0.33
##--------------------------------------------------------------------------------##
## Group means:
## Sepal.Length Sepal.Width Petal.Width Petal.Length
## setosa 5.0 3.4 0.25 1.5
## versicolor 5.9 2.8 1.33 4.3
## virginica 6.6 3.0 2.03 5.6
##--------------------------------------------------------------------------------##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.83 0.024
## Sepal.Width 1.53 2.165
## Petal.Width -2.81 2.839
## Petal.Length -2.20 -0.932
##--------------------------------------------------------------------------------##
## Proportion of trace:
## The most important part
## It shows how many discriminant axes are there and which is the most important
## in order to distinguish the taxa
## The axes are like this: a*x1 + b*x2 We get a & b from the Coefficients of
## linear discriminants
## (here: 0.83*Sepal.Length + 1.53*Sepal.Width - 2.81*Petal.Width -
## 2.20*Petal.Length for the first axis)
## We can see them with: Lec.lda$scaling
##--------------------------------------------------------------------------------##
## LD1 LD2
## 0.9912 0.0088
## LD1 is the most important as it explains 99,12%
##--------------------------------------------------------------------------------##2.1 Σημαντικότητα του κάθε διαγνωστικού άξονα σε σχέση με τους υπόλοιπους
Θέλουμε η αναλογία των τιμών αυτών, να είναι όσο το δυνατόν μεγαλύτερη, καθώς υποδηλώνουν τη σχετική σημαντικότητα κάθε άξονα. Όσο μεγαλύτερη η τιμή, τόσο πιο διαγνωστικός ο κάθε άξονας:
## [1] 48.6 4.6
2.2 Συσχέτιση των μορφολογικών γνωρισμάτων με τον κάθε διαγνωστικό άξονα
Για την εργασία σας, εάν έχετε περισσότερους από 2 διαγνωστικούς άξονες, θα χρειαστεί να τους προσθέσετε:
LD1 <- predict(Lec.lda)$x[, 1]
LD2 <- predict(Lec.lda)$x[, 2]
## We will see which variable is the most important for each axis We are
## interested in the absolute value
cor(iris[, 1:4], LD1) ## Change to mydata## [,1]
## Sepal.Length -0.79
## Sepal.Width 0.53
## Petal.Length -0.98
## Petal.Width -0.97
## [,1]
## Sepal.Length 0.218
## Sepal.Width 0.758
## Petal.Length 0.046
## Petal.Width 0.223
2.3 Πίνακας Ορθής Ταξινόμησης
Σε αυτό το σημείο, θέλουμε να δούμε πόσο καλά διακρίνουν οι μορφολογικοί χαρακτήρες οι οποίοι αναδείχθηκαν ως σημαντικοί από την ανάλυση μας, τα διάφορα taxa μεταξύ τους:
2.4 Σημαντικότητα του κάθε μορφολογικού χαρακτήρα
Τώρα μπορούμε να διαπιστώσουμε πόσο σημαντικός είναι ο εκάστοτε μορφολογικός χαρακτήρας:
## Change to mydata
## The smaller the Wilk's lambda, the better
klaR::greedy.wilks(Species ~ Sepal.Length + Sepal.Width + Petal.Width + Petal.Length,
data = iris, niveau = 0.1)## Formula containing included variables:
##
## Species ~ Petal.Length + Sepal.Width + Petal.Width + Sepal.Length
## <environment: 0x0000000015026cd0>
##
##
## Values calculated in each step of the selection procedure:
##
## vars Wilks.lambda F.statistics.overall p.value.overall
## 1 Petal.Length 0.059 1180 2.9e-91
## 2 Sepal.Width 0.037 307 2.8e-103
## 3 Petal.Width 0.025 258 5.0e-113
## 4 Sepal.Length 0.023 199 1.4e-112
## F.statistics.diff p.value.diff
## 1 1180.2 2.9e-91
## 2 43.0 2.0e-15
## 3 34.6 5.1e-13
## 4 4.7 1.0e-02
2.5 Γραφική απεικόνιση των αποτελεσμάτων
s <- predict(Lec.lda)$x
s <- as.data.frame(s)
Taxon <- factor(iris$Species) ## ## Change to mydata
s <- cbind(Taxon, s)
head(s, 1)| Taxon | LD1 | LD2 |
|---|---|---|
| Iris setosa | 8.1 | 0.3 |
## [1] "Taxon" "LD1" "LD2"
## [1] "G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/EKPA/2021"
ggsave("Iris LD1.png", width = 20, height = 20, units = "cm", dpi = 600)
ggplot(s, aes(x = LD2)) + geom_histogram(fill = "white", colour = "black") + facet_grid(Taxon ~
.)ggsave("Iris LD2.png", width = 20, height = 20, units = "cm", dpi = 600)
a <- s$LD1
b <- s$LD2
TAXON <- iris$Species ## ## Change to mydata
g <- data.frame(TAXON, a, b) ## Create a dataframe that contains all the above
## Another dataframe containing the centroids for each species and each
## discriminant axis
centroids <- aggregate(cbind(a, b) ~ TAXON, g, mean)
ggplot(g, aes(a, b, color = TAXON, shape = TAXON)) + geom_point(size = 3) + stat_ellipse(geom = "polygon",
alpha = 1/2, aes(fill = TAXON)) + xlab("LD1") + ylab("LD2") + theme(axis.title = element_text(face = "bold")) +
stat_summary(data = centroids, fun.y = "mean", geom = "point", shape = 22, size = 4,
colour = "black", aes(fill = TAXON)) +
## Shape of each species
## If you have more than 3 species you must add some values (range: 10-20)
scale_shape_manual(values = c(15, 16, 17))3 Εργασία για το σπίτι
Σας παρακαλώ θερμά να στείλετε την εργασία σας σε μορφή PDF σε αυτό το e-mail και να απαντήσετε στα ακόλουθα ερωτήματα, αφού έχετε περιγράψει το σετ δεδομένων το οποίο έχετε χρησιμοποιήσει:
1. Πόσα είδη έχει;
2. Πόσα δείγματα έχει το κάθε είδος;
3. Πόσους μορφολογικούς χαρακτήρες περιέχει το σετ δεδομένων σας;
4. Με βάση τα αποτελέσματα σας, να φτιάξετε έναν πίνακα ο οποίος να περιέχει τη μέση, την ελάχιστη και την μέγιστη τιμή για κάθε μορφολογικό χαρακτήρα κάθε είδους, μαζί με την τυπική απόκλιση.
5. Ποιες μεταβλητές δεν εμφανίζουν κανονική κατανομή;
6. Ποιες μεταβλητές εμφανίζουν μεγάλο (\(>0.8\)) συντελεστή συσχέτισης;
7. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.8\), θα τις χρησιμοποιήσετε όλες;
8. Πόσοι διαγνωστικοί άξονες εμφανίστηκαν στην ανάλυση σας; Ποιο είναι το ποσοστό της ποικιλότητας που εξηγεί ο καθένας εξ αυτών; Ποιος είναι ο σημαντικότερος και πώς μπορείτε να είστε βέβαιοι για την απάντηση σας;
9. Ποιος μορφολογικός χαρακτήρας σχετίζεται περισσότερο με κάθε εξ αυτών των αξόνων;
10. Δημιουργήστε τον πίνακα ορθής ταξινόμησης. Ποιο είναι το ποσοστό ορθής ταξινόμησης για κάθε taxon το οποίο έχετε στην ανάλυση σας; Υπάρχει ποσοστό αλληλοεπικάλυψης μεταξύ κάποιων taxa; Εάν ναι, που θεωρείτε ότι οφείλεται;
11. Ποιος μορφολογικός χαρακτήρας είναι ο σημαντικότερος για την διάκριση των taxa που έχετε στην ανάλυση σας;
12. Δημιουργήστε ένα γράφημα το οποίο να απεικονίζει γραφικά τα αποτελέσματα της ανάλυσης σας. Τι παρατηρείτε; Συμφωνούν τα αποτελέσματα του πίνακα ορθής ταξινόμησης με το εν λόγω γράφημα;
4 Appendix: R code
##===========================================================================##
## Do NOT run this code
##===========================================================================##
## Global options
library(knitr)
library(rmdformats)
## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
cache=TRUE,
prompt=FALSE,
tidy=TRUE,
# comment=NA,
message=FALSE,
warning=FALSE,
eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
pre {
max-height: 400px;
overflow-y: auto;
}
pre[class] {
max-height: 400px;
}
## CRAN repository:
install.packages("easypackages", repos = "http://cran.us.r-project.org")
library(easypackages)
libraries("rio", "car", "MASS", "psych", "leaps", "ggplot2", "klaR",
"DiscriMiner","dplyr", "car", "MASS", "corrgram", "tidyverse", "plyr")
getwd()
#setwd('E:/Academic/Lessons/Labs/Labs')
data(iris)
names(iris)
head(iris)
str(iris)
iris$Species
iris$Species <- revalue(iris$Species, c("setosa"="Iris setosa",
"virginica"="Iris virginica", "versicolor"="Iris versicolor"))
options(digits = 2)
iris %>% split(.$Species) %>% map(~summary(.))
DiscriMiner::groupStds(iris[,1:4], iris[,5])
## Violin plot
ggplot(iris, aes(Species, Petal.Length)) +
geom_violin(adjust = 4) +
geom_boxplot(width = 0.1,
fill = "black",
outlier.colour = NA) +
stat_summary(fun.y = median,
geom = "point",
fill = "white",
shape = 21,
size = 2.5) +
xlab("Taxon") +
ylab("Petal Length") +
theme(axis.title = element_text(face = "bold"),
axis.text.x = element_text(face = "italic",
colour = "black"))
getwd()
ggsave("Petal length.png",
width = 20,
height = 20,
units = "cm",
dpi = 600)
attach(iris) ## For your assignment, change that to mydata
iris <- iris[apply(iris,1,function(row) all(!is.na(row))),] ## Do the same here
set.seed(2934) ## Change that if you want to whichever number you like
cv.grp <- sample(1:10, nrow(iris),
replace = TRUE) ## Remember to change the name
cv.lda <- c() ##A function that splits our data to 2 independent subsets, one to
## train our model and one to test how well our model behaves
for(i in 1:10) {
train.ind <- which(cv.grp!=i)
test.ind <- which(cv.grp==i)
model.train <- lda(Species ~ Sepal.Length +## Change the variable names
## according to the ones that
## appear in your data
Sepal.Width + Petal.Width + Petal.Length, ## Do the same here
data = iris, subset = train.ind) ## Change to mydata
cv.lda[test.ind] <- predict(model.train,
newdata = iris[test.ind,])$class ## Do the same here
}
mean(cv.lda==as.numeric(iris$Species)) ## And here
Lec.lda <- lda(Species ~ Sepal.Length + ## This is the test point
Sepal.Width + Petal.Width + Petal.Length, ## Change the variable names
data = iris) ## Change to mydata
plot(Lec.lda,col = c("green","blue","red", "black", "orange")
[iris$Species], pch = 19) ## Change to mydata
Lec.lda
## Let's say what we did
##--------------------------------------------------------------------------------##
# Call:
# lda(Species ~ Sepal.Length + Sepal.Width + Petal.Width + Petal.Length,
# data = iris)
##--------------------------------------------------------------------------------##
## How much of the data lie in each species
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.33 0.33 0.33
##--------------------------------------------------------------------------------##
## Group means:
## Sepal.Length Sepal.Width Petal.Width Petal.Length
## setosa 5.0 3.4 0.25 1.5
## versicolor 5.9 2.8 1.33 4.3
## virginica 6.6 3.0 2.03 5.6
##--------------------------------------------------------------------------------##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.83 0.024
## Sepal.Width 1.53 2.165
##Petal.Width -2.81 2.839
##Petal.Length -2.20 -0.932
##--------------------------------------------------------------------------------##
## Proportion of trace:
## The most important part
## It shows how many discriminant axes are there and
## which is the most important in order to distinguish
## the taxa
## The axes are like this: a*x1 + b*x2
## We get a & b from the Coefficients of linear
## discriminants
##(here: 0.83*Sepal.Length + 1.53*Sepal.Width
## - 2.81*Petal.Width - 2.20*Petal.Length for the first axis)
## We can see them with: Lec.lda$scaling
##--------------------------------------------------------------------------------##
## LD1 LD2
## 0.9912 0.0088
## LD1 is the most important as it explains 99,12%
##--------------------------------------------------------------------------------##
SVD <- Lec.lda$svd
SVD
LD1 <- predict(Lec.lda)$x[,1]
LD2 <- predict(Lec.lda)$x[,2]
## We will see which variable is the most important for each axis
## We are interested in the absolute value
cor(iris[,1:4],LD1) ## Change to mydata
## Change the number of columns
cor(iris[,1:4],LD2) ## Do the same
x <- Lec.lda
y <- predict(x, iris) ## Change to mydata
detach(iris) ## Change to mydata - DON'T FORGET TO DO THAT
round(100*errormatrix(iris$Species, y$class, relative=T), 1) ## Change to mydata
## Change to mydata
## The smaller the Wilk's lambda, the better
klaR::greedy.wilks(Species ~ Sepal.Length + Sepal.Width +
Petal.Width + Petal.Length,
data = iris,
niveau = 0.1)
s <- predict(Lec.lda)$x
s <- as.data.frame(s)
Taxon <- factor(iris$Species) ## ## Change to mydata
s <- cbind(Taxon, s)
head(s,1)
names(s)
getwd()
ggplot(s, aes(x = LD1)) +
geom_histogram(fill = "white",
colour = "black") +
facet_grid(Taxon ~ .)
ggsave("Iris LD1.png", width = 20, height = 20,
units = "cm", dpi = 600)
ggplot(s, aes(x = LD2)) +
geom_histogram(fill = "white",
colour = "black") +
facet_grid(Taxon ~ .)
ggsave("Iris LD2.png",
width = 20,
height = 20,
units = "cm", dpi = 600)
a <- s$LD1
b <- s$LD2
TAXON <- iris$Species ## ## Change to mydata
g <- data.frame(TAXON, a, b) ## Create a dataframe that contains all the above
## Another dataframe containing the centroids for each species and each
## discriminant axis
centroids <- aggregate(cbind(a,b) ~ TAXON,
g,
mean)
ggplot(g,aes(a, b,
color = TAXON,
shape = TAXON)) +
geom_point(size = 3) +
stat_ellipse(geom = "polygon",
alpha = 1/2,
aes(fill = TAXON)) +
xlab("LD1") +
ylab("LD2") +
theme(axis.title = element_text(face = "bold")) +
stat_summary(data = centroids,
fun.y = "mean",
geom = "point",
shape = 22,
size = 4,
colour = "black",
aes(fill = TAXON)) +
## Shape of each species
## If you have more than 3 species you must add some values (range: 10-20)
scale_shape_manual(values = c(15, 16, 17))
ggsave("Iris ellipses.png",
width = 20,
height = 20,
units = "cm",
dpi = 600)
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##Hint: Ανατρέξτε σε προηγούμενες εργαστηριακές ασκήσεις για να δείτε πώς γίνεται ή διαβάστε το manual της εν λόγω βιβλιοθήκης: readxl↩
Δεν θα χρειαστεί να το κάνετε αυτό εσείς για την εργασία σας↩
Για την εργασία σας, για να βρείτε την τυπική απόκλιση, θα χρειαστεί να τρέξετε την ακόλουθη εντολή πρώτα:
mydata.sd <- as.data.frame(mydata)- εννοείται στην περίπτωση που έχετε ονομάσει το αρχείο σας mydata↩Για να κάνετε αυτό το βήμα για την εργασία σας, θα πρέπει να αλλάξετε το όνομα του αρχείου, καθώς και τα ονόματα των μεταβλητών - Εάν έχετε απορίες, τώρα είναι η στιγμή να τις εκφράσετε↩